http://geophysique.be/2014/02/25/shaded-relief-map-in-python/
Author: Thomas Lecocq
Some standard and some less-standard imports:
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
from cartopy.io.img_tiles import GoogleTiles
from cartopy.io.srtm import srtm_composite
from osgeo import gdal
from osgeo import gdal_array
Reasons for each packages:
I am currently in Sudelfeld for a workshop organised by the colleagues of Munich University. Let's plot a 1-degree square image comprising the village we are in. The Google Maps tiles are downloaded automatically from the webservice (so an Internet connection is needed):
# Specify a region of interest, in this case, Sudelfeld Ski Resort (Germany)
lat = 47 + 40 / 60.0 + 30 / 3600.
lon = 12 + 3 / 60.0 + 2 / 3600.
plt.figure(figsize=(10, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
ax.set_extent([12.0, 13.0, 47.0, 48.0])
gg_tiles = GoogleTiles()
ax.add_image(gg_tiles, 10)
plt.scatter(lon, lat, marker=(5, 1), color='red', s=200)
plt.title("Welcome to Sudelfeld")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False
SRTM data can be fetched automatically like the tiles, and plotted using imshow:
plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
elev, crs, extent = srtm_composite(12, 47, 1, 1)
plt.imshow(elev, extent=extent, transform=crs,
cmap='gist_earth', origin='lower')
cb = plt.colorbar(orientation='vertical')
cb.set_label('Altitude')
plt.title("SRTM Map")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False
As visible above, there are bad data (Gaps) in the SRTM tile we downloaded. This kind of lacking data is particularly annoying, so, let's mask all values below 0:
plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
elev, crs, extent = srtm_composite(12, 47, 1, 1)
elev = np.ma.masked_less_equal(elev,0,copy=False)
plt.imshow(elev, extent=extent, transform=crs,
cmap='gist_earth', origin='lower')
cb = plt.colorbar(orientation='vertical')
cb.set_label('Altitude')
plt.title("SRTM Map")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False
This is way better... but gaps/masked data are going to be quite annoying for the slope/aspect calculations, so let's use GDAL's methods to interpolate data from around the gaps:
plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
elev, crs, extent = srtm_composite(12, 47, 1, 1)
src_ds = gdal_array.OpenArray(elev)
srcband = src_ds.GetRasterBand(1)
dstband = srcband
maskband = srcband
smoothing_iterations = 0
options = []
max_distance = 200
result = gdal.FillNodata(dstband, maskband,
max_distance, smoothing_iterations, options,
callback=None)
elev = dstband.ReadAsArray()
plt.imshow(elev, extent=extent, transform=crs,
cmap='gist_earth', origin='lower')
cb = plt.colorbar(orientation='vertical')
cb.set_label('Altitude')
plt.title("SRTM Map")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False
Now, that's way better !
A shaded relief map is computed using the slope and the aspect of the topography. We want to have a look at the map, in the morning (Sun on the East, at 30° above the horizon):
## Subplot 23: Shaded Relief Map
plt.figure(figsize=(10, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
x, y = np.gradient(elev)
slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
# -x here because of pixel orders in the SRTM tile
aspect = np.arctan2(-x, y)
altitude = np.pi / 6.
azimuth = np.pi /2.
shaded = np.sin(altitude) * np.sin(slope)\
+ np.cos(altitude) * np.cos(slope)\
* np.cos((azimuth - np.pi/2.) - aspect)
plt.imshow(shaded, extent=extent, transform=crs,
cmap='Greys', origin='lower')
plt.title("Shaded Relief Map")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False
Done !